1 Setup

1.1 Libs

library(car) 
## Loading required package: carData
library("Rmisc")
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.3.1
## Loading required package: plyr
## Warning: package 'plyr' was built under R version 4.3.1
library("ggpubr")
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.3.1
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
## 
##     mutate
library("cowplot")
## Warning: package 'cowplot' was built under R version 4.3.1
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
library("lme4")
## Warning: package 'lme4' was built under R version 4.3.1
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.3.1
## Warning in check_dep_version(): ABI version mismatch: 
## lme4 was built with Matrix ABI version 1
## Current Matrix ABI version is 0
## Please re-install lme4 from source or restore original 'Matrix' package
library("phyloseq")
library("glmmTMB")
## Warning: package 'glmmTMB' was built under R version 4.3.1
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.6.0
## Current Matrix version is 1.6.1.1
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
library("emmeans")
## Warning: package 'emmeans' was built under R version 4.3.1
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
#install.packages("multcompView")
library(multcompView)
library("bbmle")
## Loading required package: stats4
library("DHARMa")
## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library("nlme")
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
## 
##     lmList
library("bestNormalize")
## 
## Attaching package: 'bestNormalize'
## The following object is masked from 'package:MASS':
## 
##     boxcox

1.2 Data

setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Mosquito_business/Mosquito_microbes/04.microbiome_analysis/04a.diversity")

ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.less.rds")
ps.clean 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 161 taxa and 458 samples ]
## sample_data() Sample Data:       [ 458 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 161 taxa by 11 taxonomic ranks ]
#ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.rare9200.rds")
#ps.clean

samdf <- data.frame(ps.clean@sam_data)

##made below
#df.div <- read.csv("div.trim.csv",row.names=1)

2 Alpha diversity calculations (first run through)

##option to run rarefied stuff:
#ps.clean <- ps.rare

df.all <- data.frame(estimate_richness(ps.clean, split=TRUE, measures=c("Shannon","InvSimpson","Observed")))

df.all$sample_name <- rownames(df.all)

df.all.div <- merge(df.all,samdf,by="sample_name") #add sample data

#shannon diversity divided by species richness
df.all.div$even <- df.all.div$Shannon/(log(df.all.div$Observed))

df.div <- df.all.div

#adding most recent sample counts
sums <- data.frame(sample_sums(ps.clean))
colnames(sums) <- c("lib_size_trim")
sums$sample_name <- row.names(sums)

df.div <- merge(df.div,sums,by="sample_name")

#write.csv(df.div,"div.trim.csv")

3 Analyzing alpha metrics

3.1 All data

Very strange/interesting pattern where the water types that shouldn’t have a lot of diversity have a lot of OTUs, but not a lot of reads…

ggplot(df.div,aes(x=type,y=Observed,color=infusion))+
  geom_boxplot()+
  geom_jitter()

ggplot(df.div,aes(x=type,y=lib_size_trim,color=infusion))+
  geom_boxplot()+
  geom_jitter(position=position_jitterdodge())

3.2 Plotz

#just experimental types
df.div.exp <- subset(df.div,type=="Microbial Water"|type=="A.albopictus")

df.div.exp$inf.temp <- paste0(df.div.exp$infusion,df.div.exp$temperature)

df.div.exp$type <- sub("A.albopictus","Mosquitoes",df.div.exp$type)
df.div.exp$type <- sub("Microbial Water","Meso. water",df.div.exp$type)

df.div.exp$type <- factor(df.div.exp$type,levels=c("Mosquitoes","Meso. water"))

3.2.1 Richness

df.div.exp.se <- summarySE(df.div.exp,measurevar="Observed",groupvars=c("inf.temp","infusion","temperature","type"))

##dots & error bars
gg.rich <- ggplot(df.div.exp,aes(x=infusion,y=Observed,fill=inf.temp,shape=inf.temp))+
  geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  ylim(4,34)+
  facet_wrap(~type)+
  theme_cowplot()+
  ylab("ASV richness")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))
gg.rich

##water by time
df.div.water <- subset(df.div.exp,type=="Meso. water")
gg.mw.rich.time <- ggplot(df.div.water,aes(x=infusion,y=Observed,fill=inf.temp))+
  geom_boxplot()+
  #geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  #geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  #geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  facet_wrap(~day)+
  theme_cowplot()+
  ylab("ASV richness")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  ggtitle("Mesocosm water")
gg.mw.rich.time

3.2.2 Simpson’s

df.div.exp.se <- summarySE(df.div.exp,measurevar="InvSimpson",groupvars=c("inf.temp","infusion","temperature","type"))

##dots & error bars
gg.simp <- ggplot(df.div.exp,aes(x=infusion,y=InvSimpson,fill=inf.temp,shape=inf.temp))+
  geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  geom_errorbar(data=df.div.exp.se,aes(ymax=InvSimpson+se,ymin=InvSimpson-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  facet_wrap(~type)+
  theme_cowplot()+
  ylab("Simpson's (inv.)")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))
gg.simp

gg.mw.simp.time <- ggplot(df.div.water,aes(x=infusion,y=InvSimpson,fill=inf.temp))+
  geom_boxplot()+
  #geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  #geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  #geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  facet_wrap(~day)+
  theme_cowplot()+
  ylab("Simpson's index (inv.)")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  ggtitle("Mesocosm water")

3.2.3 Multi-panel (Figure 2)

gg.div <- ggarrange(gg.rich,gg.simp,common.legend=T,legend="right",labels=c("(a)","(b)"))
gg.div

##just water things
ggarrange(gg.mw.rich.time,gg.mw.simp.time,common.legend=T,legend="right",nrow=2)

#ggsave("div.water.time.png",width=6,height=6)

3.3 Statz

3.3.1 Richness

3.3.1.1 Water richness stats

df.div.mw <- subset(df.div.exp,type=="Meso. water")

df.div.mw$day <- as.factor(df.div.mw$day)
str(df.div.mw)
## 'data.frame':    211 obs. of  30 variables:
##  $ sample_name      : chr  "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
##  $ Observed         : num  25 28 22 28 31 31 28 29 27 21 ...
##  $ Shannon          : num  1.63 1.53 1.58 1.53 1.75 ...
##  $ InvSimpson       : num  2.67 2.4 2.44 2.32 2.93 ...
##  $ orgname          : chr  "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
##  $ mesocosm         : chr  "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
##  $ type             : Factor w/ 2 levels "Mosquitoes","Meso. water": 2 2 2 2 2 2 2 2 2 2 ...
##  $ stage            : chr  "Water" "Water" "Water" "Water" ...
##  $ sex              : chr  "" "" "" "" ...
##  $ date.collected   : chr  "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
##  $ day              : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
##  $ abdomen.length.mm: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ wing.length      : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ hindleg.length   : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ special          : chr  "" "" "" "" ...
##  $ dispersal        : chr  "N" "N" "D" "D" ...
##  $ temperature      : chr  "C" "C" "C" "C" ...
##  $ infusion         : chr  "OL" "OL" "OL" "OL" ...
##  $ raw_miseq_reads  : int  30391 39927 18829 58206 54860 39963 37857 45572 36841 27161 ...
##  $ X16scq           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Wolb_ct1         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Wolb_ct2         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Act_ct           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ notes            : chr  "" "" "" "" ...
##  $ newday           : chr  "early" "early" "early" "early" ...
##  $ lib_size         : num  25797 34139 15861 49298 44559 ...
##  $ is.neg           : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ even             : num  0.506 0.458 0.511 0.46 0.51 ...
##  $ lib_size_trim    : num  25797 34134 15861 49290 44524 ...
##  $ inf.temp         : chr  "OLC" "OLC" "OLC" "OLC" ...
#hist(df.div.mw$Observed)
#shapiro.test(log(df.div.mw$Observed))
shapiro.test(df.div.mw$Observed)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mw$Observed
## W = 0.9775, p-value = 0.001855
#wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mw)
#wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)

#AICtab(wrich1,wrich2,wrich3)

##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus temperature 
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus time 
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)

anova(wrich1,wrich1nd) #0.3931
## Data: df.div.mw
## Models:
## wrich1nd: Observed ~ day + temperature + infusion + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nd:     (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1nd  8 1035.4 1062.2 -509.69   1019.4                         
## wrich1    9 1036.7 1066.8 -509.33   1018.7 0.7293      1     0.3931
anova(wrich1,wrich1ni) #2.2e-16***
## Data: df.div.mw
## Models:
## wrich1ni: Observed ~ day + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1ni:     (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wrich1ni  7 1118.0 1141.5 -552.02   1104.0                             
## wrich1    9 1036.7 1066.8 -509.33   1018.7 85.373      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wrich1,wrich1nt) #0.4119
## Data: df.div.mw
## Models:
## wrich1nt: Observed ~ day + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nt:     (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1nt  8 1035.3 1062.2 -509.67   1019.3                         
## wrich1    9 1036.7 1066.8 -509.33   1018.7 0.6732      1     0.4119
anova(wrich1,wrich1nnd) #0.001874**
## Data: df.div.mw
## Models:
## wrich1nnd: Observed ~ temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nnd:     (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## wrich1nnd  7 1045.2 1068.7 -515.61   1031.2                            
## wrich1     9 1036.7 1066.8 -509.33   1018.7 12.559      2   0.001874 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(wrich1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                Chisq Df Pr(>Chisq)    
## day          12.9404  2   0.001549 ** 
## temperature   0.6742  1   0.411577    
## infusion    151.4172  2  < 2.2e-16 ***
## dispersal     0.7306  1   0.392704    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Response: Observed
# Response: Observed
#                Chisq Df Pr(>Chisq)    
# day          12.9404  2   0.001549 ** 
# temperature   0.6742  1   0.411577    
# infusion    151.4172  2  < 2.2e-16 ***
# dispersal     0.7306  1   0.392704    
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... almost tied with no dispersal, then full
##           dAIC df
## wrich1nt   0.0 8 
## wrich1nd   0.1 8 
## wrich1     1.3 9 
## wrich1nnd  9.9 7 
## wrich1ni  82.7 7
##interactions?
wrich1.int <- glmmTMB(Observed~day*infusion+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int1 <- glmmTMB(Observed~day+infusion*temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int2 <- glmmTMB(Observed~day+infusion+temperature*dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int3 <- glmmTMB(Observed~day*temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int4 <- glmmTMB(Observed~day*dispersal+temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int5 <- glmmTMB(Observed~day+dispersal*infusion+temperature+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)

anova(wrich1,wrich1.int) #2.86e-10***
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int:     (1 | mesocosm), zi=~0, disp=~1
##            Df     AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wrich1      9 1036.66 1066.8 -509.33  1018.66                             
## wrich1.int 13  994.79 1038.4 -484.40   968.79 49.863      4  3.857e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wrich1,wrich1.int1) #ns
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int1: Observed ~ day + infusion * temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int1:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1036.7 1066.8 -509.33   1018.7                         
## wrich1.int1 11 1037.2 1074.1 -507.62   1015.2 3.4154      2     0.1813
anova(wrich1,wrich1.int2) #ns
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int2: Observed ~ day + infusion + temperature * dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int2:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1036.7 1066.8 -509.33   1018.7                         
## wrich1.int2 10 1037.1 1070.7 -508.57   1017.1 1.5253      1     0.2168
anova(wrich1,wrich1.int3) #ns
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int3: Observed ~ day * temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int3:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1036.7 1066.8 -509.33   1018.7                         
## wrich1.int3 11 1039.0 1075.9 -508.52   1017.0 1.6226      2     0.4443
anova(wrich1,wrich1.int4) #ns
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int4: Observed ~ day * dispersal + temperature + infusion + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int4:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1036.7 1066.8 -509.33   1018.7                         
## wrich1.int4 11 1039.3 1076.2 -508.66   1017.3 1.3341      2     0.5132
anova(wrich1,wrich1.int5) #ns
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int5: Observed ~ day + dispersal * infusion + temperature + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int5:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1036.7 1066.8 -509.33   1018.7                         
## wrich1.int5 11 1037.8 1074.7 -507.90   1015.8 2.8599      2     0.2393
Anova(wrich1.int)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                 Chisq Df Pr(>Chisq)    
## day           17.4811  2    0.00016 ***
## infusion     167.8537  2  < 2.2e-16 ***
## temperature    0.7595  1    0.38348    
## dispersal      0.7627  1    0.38248    
## day:infusion  58.9948  4  4.717e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Observed
#                 Chisq Df Pr(>Chisq)    
# day           17.4811  2    0.00016 ***
# infusion     167.8537  2  < 2.2e-16 ***
# temperature    0.7595  1    0.38348    
# dispersal      0.7627  1    0.38248    
# day:infusion  58.9948  4  4.717e-12 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#tapply(df.div.mw$log(Observed), df.div.mw$infusion, mean)
#tapply(df.div.mw$log(Observed), df.div.mw$day, mean)

wrich1.red <- glmmTMB(Observed~day+infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.red.int <- glmmTMB(Observed~day*infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)

anova(wrich1.red,wrich1.red.int) #sig 4e-10***
## Data: df.div.mw
## Models:
## wrich1.red: Observed ~ day + infusion + offset(log(lib_size_trim)) + (1 | , zi=~0, disp=~1
## wrich1.red:     mesocosm), zi=~0, disp=~1
## wrich1.red.int: Observed ~ day * infusion + offset(log(lib_size_trim)) + (1 | , zi=~0, disp=~1
## wrich1.red.int:     mesocosm), zi=~0, disp=~1
##                Df     AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wrich1.red      7 1034.08 1057.5 -510.04  1020.08                             
## wrich1.red.int 11  992.32 1029.2 -485.16   970.32 49.758      4  4.056e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##model checking
shapiro.test(residuals(wrich1.red.int)) #awesome
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(wrich1.red.int)
## W = 0.9921, p-value = 0.3128
qqnorm(residuals(wrich1.red.int))
qqline(residuals(wrich1.red.int), col="red")

#plotResiduals(wrich1.int)
#wrich1.red.int.resid <- simulateResiduals(fittedModel = wrich1.red.int, plot = T)

##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.int.em <- emmeans(wrich1.int,~infusion*temperature)

multcomp::cld(wrich1.int.em)
##  infusion temperature emmean    SE  df lower.CL upper.CL .group
##  SW       C             19.1 0.351 198     18.5     19.8  1    
##  SW       H             19.5 0.353 198     18.8     20.1  1    
##  SG       C             21.2 0.355 198     20.5     21.9   2   
##  SG       H             21.5 0.357 198     20.8     22.2   2   
##  OL       C             24.7 0.351 198     24.0     25.4    3  
##  OL       H             25.0 0.353 198     24.3     25.7    3  
## 
## Results are averaged over the levels of: day, dispersal 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
 # infusion temperature emmean    SE  df lower.CL upper.CL .group
 # SW       C             19.1 0.351 198     18.5     19.8  1    
 # SW       H             19.5 0.353 198     18.8     20.1  1    
 # SG       C             21.2 0.355 198     20.5     21.9   2   
 # SG       H             21.5 0.357 198     20.8     22.2   2   
 # OL       C             24.7 0.351 198     24.0     25.4    3  
 # OL       H             25.0 0.353 198     24.3     25.7    3  

3.3.1.2 Water richness stats - rarefied

Same as above but without offset for library size. Not knitting this part, just checking out in R

df.div.mw <- subset(df.div.exp,type=="Meso. water")

df.div.mw$day <- as.factor(df.div.mw$day)
str(df.div.mw)

#hist(df.div.mw$Observed)
#shapiro.test(log(df.div.mw$Observed))
shapiro.test(df.div.mw$Observed)

#wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal)+(1|mesocosm),family="poisson",data=df.div.mw)
#wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal)+(1|mesocosm),family="compois",data=df.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)

#AICtab(wrich1,wrich2,wrich3)

##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+(1|mesocosm), data=df.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+(1|mesocosm), data=df.div.mw)
##minus temperature 
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+(1|mesocosm), data=df.div.mw)
##minus time 
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+(1|mesocosm), data=df.div.mw)

anova(wrich1,wrich1nd) #0.466
anova(wrich1,wrich1ni) #2.2e-16***
anova(wrich1,wrich1nt) #0.4789
anova(wrich1,wrich1nnd) #0.006693**

Anova(wrich1)
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Observed
#                Chisq Df Pr(>Chisq)    
# day          10.2547  2   0.005932 ** 
# temperature   0.5020  1   0.478604    
# infusion    133.7700  2  < 2.2e-16 ***
# dispersal     0.5322  1   0.465695    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... almost tied with no dispersal, then full

##interactions?
wrich1.int <- glmmTMB(Observed~day*infusion+temperature+dispersal+(1|mesocosm),data=df.div.mw)
wrich1.int1 <- glmmTMB(Observed~day+infusion*temperature+dispersal+(1|mesocosm),data=df.div.mw)
wrich1.int2 <- glmmTMB(Observed~day+infusion+temperature*dispersal+(1|mesocosm),data=df.div.mw)
wrich1.int3 <- glmmTMB(Observed~day*temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
wrich1.int4 <- glmmTMB(Observed~day*dispersal+temperature+infusion+(1|mesocosm),data=df.div.mw)
wrich1.int5 <- glmmTMB(Observed~day+dispersal*infusion+temperature+(1|mesocosm),data=df.div.mw)

anova(wrich1,wrich1.int) #4.35e-11***
anova(wrich1,wrich1.int1) #ns
anova(wrich1,wrich1.int2) #ns
anova(wrich1,wrich1.int3) #ns
anova(wrich1,wrich1.int4) #ns
anova(wrich1,wrich1.int5) #ns

Anova(wrich1.int)
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Observed
#                 Chisq Df Pr(>Chisq)    
# day           13.9847  2  0.0009189 ***
# infusion     154.5860  2  < 2.2e-16 ***
# temperature    0.5795  1  0.4464988    
# dispersal      0.5809  1  0.4459484    
# day:infusion  64.8562  4  2.759e-13 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#tapply(df.div.mw$log(Observed), df.div.mw$infusion, mean)
#tapply(df.div.mw$log(Observed), df.div.mw$day, mean)

wrich1.red <- glmmTMB(Observed~day+infusion+(1|mesocosm),data=df.div.mw)
wrich1.red.int <- glmmTMB(Observed~day*infusion+(1|mesocosm),data=df.div.mw)

anova(wrich1.red,wrich1.red.int) #sig 4e-11**

##model checking
shapiro.test(residuals(wrich1.red.int)) #awesome
qqnorm(residuals(wrich1.red.int))
qqline(residuals(wrich1.red.int), col="red")

#plotResiduals(wrich1.int)
#wrich1.red.int.resid <- simulateResiduals(fittedModel = wrich1.red.int, plot = T)

##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.int.em <- emmeans(wrich1.int,~infusion*temperature)

multcomp::cld(wrich1.int.em)

3.3.1.3 Mosquito richness stats

df.div.mq <- subset(df.div.exp,type=="Mosquitoes")

hist(df.div.mq$Observed)

shapiro.test(df.div.mq$Observed)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mq$Observed
## W = 0.89464, p-value = 1.697e-10
##different stat families
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
mrich2 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mq)
mrich3 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="genpois",data=df.div.mq)
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mq)
#mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="nbinom1",data=df.div.mq)
#mrich5 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="nbinom2",data=df.div.mq)

AICtab(mrich1,mrich2,mrich3,mrich4) #gaussian again
##        dAIC df
## mrich1  0.0 10
## mrich2 47.1 10
## mrich4 48.3 9 
## mrich3 48.3 10
##full model
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##no dispersal 
mrich1nd <- glmmTMB(Observed~newday+temperature+infusion+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus infusion 
mrich1ni <- glmmTMB(Observed~newday+temperature+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus temperature 
mrich1nt <- glmmTMB(Observed~newday+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus time 
mrich1nnd <- glmmTMB(Observed~temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mq)
##minus sex
mrich1ns <- glmmTMB(Observed~newday+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)

Anova(mrich1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##              Chisq Df Pr(>Chisq)
## newday      1.8273  2     0.4011
## temperature 0.2244  1     0.6357
## infusion    0.4926  2     0.7817
## dispersal   1.4645  1     0.2262
## sex         0.3992  1     0.5275
shapiro.test(residuals(mrich1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mrich1)
## W = 0.90866, p-value = 1.328e-09
qqnorm(residuals(mrich1))
qqline(residuals(mrich1), col="red")

mrich1.resid <- simulateResiduals(fittedModel = mrich1, plot = T) 

##posthoc things
mrich1.em <- emmeans(mrich1,~infusion*temperature)
multcomp::cld(mrich1.em)
##  infusion temperature emmean    SE  df lower.CL upper.CL .group
##  OL       H             8.34 0.417 185     7.51     9.16  1    
##  SG       H             8.49 0.447 185     7.61     9.37  1    
##  OL       C             8.56 0.340 185     7.89     9.23  1    
##  SG       C             8.72 0.490 185     7.75     9.68  1    
##  SW       H             8.77 0.479 185     7.82     9.71  1    
##  SW       C             8.99 0.564 185     7.88    10.10  1    
## 
## Results are averaged over the levels of: newday, dispersal, sex 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
##no differences

3.3.1.4 Mosquito richness stats - rarefied

Same as above but without offset for library size

df.div.mq <- subset(df.div.exp,type=="Mosquitoes")

hist(df.div.mq$Observed)
shapiro.test(df.div.mq$Observed)

##different stat families
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),data=df.div.mq)
mrich2 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
mrich3 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="genpois",data=df.div.mq)
mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="poisson",data=df.div.mq)
#mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="nbinom1",data=df.div.mq)
#mrich5 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="nbinom2",data=df.div.mq)

AICtab(mrich1,mrich2,mrich3,mrich4) #compois 

##full model
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##no dispersal 
mrich1nd <- glmmTMB(Observed~newday+temperature+infusion+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus infusion 
mrich1ni <- glmmTMB(Observed~newday+temperature+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus temperature 
mrich1nt <- glmmTMB(Observed~newday+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus time 
mrich1nnd <- glmmTMB(Observed~temperature+infusion+dispersal+sex+(1|mesocosm),family="compois", data=df.div.mq)
##minus sex
mrich1ns <- glmmTMB(Observed~newday+temperature+infusion+dispersal+(1|mesocosm),family="compois",data=df.div.mq)

Anova(mrich1)

shapiro.test(residuals(mrich1))
qqnorm(residuals(mrich1))
qqline(residuals(mrich1), col="red")

mrich1.resid <- simulateResiduals(fittedModel = mrich1, plot = T) 

##posthoc things
mrich1.em <- emmeans(mrich1,~infusion*temperature)
multcomp::cld(mrich1.em)
##no differences

3.3.2 Simpson’s

3.3.2.1 Water Simpson’s stats

hist(df.div.mw$InvSimpson)

shapiro.test(df.div.mw$InvSimpson)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mw$InvSimpson
## W = 0.96756, p-value = 9.003e-05
hist(log(df.div.mw$InvSimpson))

shapiro.test(log(df.div.mw$InvSimpson))
## 
##  Shapiro-Wilk normality test
## 
## data:  log(df.div.mw$InvSimpson)
## W = 0.95673, p-value = 5.166e-06
##normalizing because residuals bad below
bestNormalize(df.div.mw$InvSimpson) #orderNorm
## Best Normalizing transformation with 211 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.5685
##  - Box-Cox: 1.4779
##  - Center+scale: 1.4303
##  - Double Reversed Log_b(x+a): 1.5946
##  - Exp(x): 3.4885
##  - Log_b(x+a): 1.6482
##  - orderNorm (ORQ): 1.2276
##  - sqrt(x + a): 1.4479
##  - Yeo-Johnson: 1.4679
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 211 nonmissing obs and no ties 
##  - Original quantiles:
##    0%   25%   50%   75%  100% 
## 1.359 1.947 2.669 3.225 5.130
df.div.mw$simp.norm <- bestNormalize(df.div.mw$InvSimpson)$x.t
shapiro.test(df.div.mw$simp.norm)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mw$simp.norm
## W = 0.99978, p-value = 1
##just replaced InvSimpson with simp.norm
##full model
wsimp1 <- glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##no dispersal
wsimp1nd <- glmmTMB(simp.norm~day+temperature+infusion+(1|mesocosm),data=df.div.mw)
##minus infusion
wsimp1ni <- glmmTMB(simp.norm~day+temperature+dispersal+(1|mesocosm),data=df.div.mw)
##minus temperature 
wsimp1nt <- glmmTMB(simp.norm~day+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##minus time 
wsimp1nnd <- glmmTMB(simp.norm~temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)

Anova(wsimp1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##                Chisq Df Pr(>Chisq)    
## day          49.7901  2  1.542e-11 ***
## temperature  41.3001  1  1.306e-10 ***
## infusion    229.1055  2  < 2.2e-16 ***
## dispersal     0.0091  1      0.924    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: simp.norm
#                Chisq Df Pr(>Chisq)    
# day          49.7901  2  1.542e-11 ***
# temperature  41.3001  1  1.306e-10 ***
# infusion    229.1055  2  < 2.2e-16 ***
# dispersal     0.0091  1      0.924    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##~same results with normalized results
##same results with rarefied

AICtab(wsimp1,wsimp1nd,wsimp1ni,wsimp1nt,wsimp1nnd)
##           dAIC  df
## wsimp1nd    0.0 8 
## wsimp1      2.0 9 
## wsimp1nt   32.6 8 
## wsimp1nnd  40.7 7 
## wsimp1ni  101.8 7
##no dispersal good, then full model

##interaction?
wsimp1.int1 <- glmmTMB(simp.norm~day*infusion*temperature+(1|mesocosm), data=df.div.mw)
wsimp1.int2 <- glmmTMB(simp.norm~day*infusion+temperature+(1|mesocosm), data=df.div.mw) #best AIC
wsimp1.int3 <- glmmTMB(simp.norm~day+infusion*temperature+(1|mesocosm), data=df.div.mw)
wsimp1.int4 <- glmmTMB(simp.norm~day*temperature+infusion+(1|mesocosm), data=df.div.mw)

AICtab(wsimp1.int1,wsimp1.int2,wsimp1.int3,wsimp1.int4,wsimp1nd) #int better
##             dAIC df
## wsimp1.int2  0.0 12
## wsimp1.int1  5.6 20
## wsimp1.int3 40.5 10
## wsimp1nd    40.6 8 
## wsimp1.int4 44.3 10
#interaction with day & infusion

Anova(wsimp1.int2) #all 0.001 sig level, same with rare
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##                Chisq Df Pr(>Chisq)    
## day           70.733  2  4.370e-16 ***
## infusion     227.491  2  < 2.2e-16 ***
## temperature   39.606  1  3.107e-10 ***
## day:infusion  58.541  4  5.874e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Response: simp.norm
#                Chisq Df Pr(>Chisq)    
# day           70.733  2  4.370e-16 ***
# infusion     227.491  2  < 2.2e-16 ***
# temperature   39.606  1  3.107e-10 ***
# day:infusion  58.541  4  5.874e-12 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

##all of the below was bad before normalizing
shapiro.test(residuals(wsimp1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(wsimp1)
## W = 0.96148, p-value = 1.724e-05
qqnorm(residuals(wsimp1))
qqline(residuals(wsimp1), col="red")

wsimp1.resid <- simulateResiduals(fittedModel = wsimp1, plot = T) 

##posthoc things
wsimp1.em <- emmeans(wsimp1,~infusion*temperature)
multcomp::cld(wsimp1.em)
##  infusion temperature emmean     SE  df lower.CL upper.CL .group
##  SG       C           -1.255 0.0934 202   -1.439  -1.0708  1    
##  SG       H           -0.656 0.0939 202   -0.841  -0.4712   2   
##  OL       C           -0.104 0.0927 202   -0.287   0.0786    3  
##  SW       C            0.442 0.0927 202    0.259   0.6244     4 
##  OL       H            0.495 0.0932 202    0.311   0.6783     4 
##  SW       H            1.040 0.0932 202    0.857   1.2240      5
## 
## Results are averaged over the levels of: day, dispersal 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
 # infusion temperature emmean     SE  df lower.CL upper.CL .group
 # SG       C           -1.255 0.0934 202   -1.439  -1.0708  1    
 # SG       H           -0.656 0.0939 202   -0.841  -0.4712   2   
 # OL       C           -0.104 0.0927 202   -0.287   0.0786    3  
 # SW       C            0.442 0.0927 202    0.259   0.6244     4 
 # OL       H            0.495 0.0932 202    0.311   0.6783     4 
 # SW       H            1.040 0.0932 202    0.857   1.2240      5

3.3.2.2 Mosquito Simpson’s

From manuscript pre-revisions: “In contrast, the Simpson’s index, which integrates evenness into a measure of alpha diversity, of the adult mosquito microbiome was lower in male mosquitoes relative to female mosquitoes (p<0.0001). This trend was associated with the increased dominance of wAlbB in males. The Simpson’s index of the adult mosquito microbiome did not vary significantly with the dispersal, aquatic chemistry, temperature, or time of emergence.”

hist(df.div.mq$InvSimpson)

df.div.mq$simp.norm <- bestNormalize(df.div.mq$InvSimpson)$x.t #orderNorm
hist(df.div.mq$simp.norm)

##full model
msimp.all <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without sex
msimp.nosex <- glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm), data=df.div.mq)
##without dispersal
msimp.nodis <- glmmTMB(simp.norm~newday+temperature+infusion+sex+(1|mesocosm), data=df.div.mq)
##without infusion
msimp.noinf <- glmmTMB(simp.norm~newday+temperature+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without temperature
msimp.notem <- glmmTMB(simp.norm~newday+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without time
msimp.notim <- glmmTMB(simp.norm~temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)

Anova(msimp.all)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##               Chisq Df Pr(>Chisq)    
## newday       2.3823  2    0.30387    
## temperature  0.3062  1    0.58000    
## infusion     2.7438  2    0.25363    
## dispersal    2.8697  1    0.09026 .  
## sex         67.4477  1    < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: simp.norm
#               Chisq Df Pr(>Chisq)    
# newday       2.2911  2    0.31805    
# temperature  0.2987  1    0.58469    
# infusion     2.7011  2    0.25909    
# dispersal    2.8332  1    0.09234 .  
# sex         67.5319  1    < 2e-16 ***

AICtab(msimp.all,msimp.nosex,msimp.nodis,msimp.noinf,msimp.notem,msimp.notim)
##             dAIC df
## msimp.notem  0.0 9 
## msimp.notim  0.1 8 
## msimp.noinf  0.4 8 
## msimp.all    1.7 10
## msimp.nodis  2.4 9 
## msimp.nosex 59.0 8
#no time & no temp tied for best, then no infusion

##all of the below was bad before normalizing
shapiro.test(residuals(msimp.all))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(msimp.all)
## W = 0.98791, p-value = 0.0961
qqnorm(residuals(msimp.all))
qqline(residuals(msimp.all), col="red")

msimp.all.resid <- simulateResiduals(fittedModel = msimp.all, plot = T) 

##posthoc things
msimp.all.em <- emmeans(msimp.all,~infusion*temperature)
multcomp::cld(msimp.all.em)
##  infusion temperature  emmean    SE  df lower.CL upper.CL .group
##  SG       H           -0.2061 0.152 185   -0.506   0.0939  1    
##  SG       C           -0.1150 0.167 185   -0.445   0.2150  1    
##  OL       H           -0.0782 0.143 185   -0.361   0.2045  1    
##  OL       C            0.0129 0.115 185   -0.215   0.2406  1    
##  SW       H            0.1484 0.162 185   -0.172   0.4682  1    
##  SW       C            0.2394 0.193 185   -0.141   0.6197  1    
## 
## Results are averaged over the levels of: newday, dispersal, sex 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
#no diffs

3.4 Infusion water things

df.div.inf <- subset(df.div,type=="Infusion water")

df.div.inf.se <- summarySE(df.div.inf,measurevar="Observed",groupvars=c("infusion"))
df.div.inf.se
##   infusion N  Observed       sd        se        ci
## 1       OL 3 45.666667 3.214550 1.8559215  7.985386
## 2       SG 3 23.666667 6.350853 3.6666667 15.776393
## 3       SW 3  4.333333 1.154701 0.6666667  2.868435
# gg.iw.rich <- ggplot(df.div.inf,aes(x=infusion,y=Observed,fill=infusion))+
#   geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
#   geom_errorbar(data=df.div.inf.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
#   geom_point(data=df.div.inf.se,size=2.5,position=position_dodge(width=0.6))+
#   #geom_boxplot(outlier.shape=NA,alpha=0.5)+
#   theme_cowplot()+
#   ggtitle("Infusion water")
# gg.iw.rich

4 Diversity correlations

4.1 Diversity x day

Doesn’t look interesting

ggplot(df.div.mq,aes(x=day,y=Observed))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(df.div.mq,aes(x=day,y=InvSimpson))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

4.2 Diversity - mesocosm vs. mosquito

library(reshape)
## 
## Attaching package: 'reshape'
## The following object is masked from 'package:Matrix':
## 
##     expand
## The following object is masked from 'package:cowplot':
## 
##     stamp
## The following objects are masked from 'package:plyr':
## 
##     rename, round_any
meso.corr <- merge(df.div.mq,df.div.mw,by="mesocosm")
ggplot(meso.corr,aes(x=InvSimpson.x,y=InvSimpson.y))+
  geom_point()+
  #facet_wrap(~day.y)+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(meso.corr,aes(x=Observed.x,y=Observed.y))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'